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We calculate the gravitational self force acting on a pointlike particle of mass /i, set in a circular 
geodesic orbit around a Schwarzschild black hole. Our calculation is done in the Lorenz gauge: For 
given orbital radius, we first solve directly for the Lorenz-gauge metric perturbation using numerical 
evolution in the time domain; We then compute the (finite) back-reaction force from each of the 
multipole modes of the perturbation; Finally, we apply the "mode sum" method to obtain the total, 
physical self force. The temporal component of the self force (which is gauge invariant) describes 
the dissipation of orbital energy through gravitational radiation. Our results for this component are 
consistent, to within the computational accuracy, with the total flux of gravitational-wave energy 
radiated to infinity and through the event horizon. The radial component of the self force (which is 
gauge dependent) is calculated here for the first time. It describes a conservative shift in the orbital 
parameters away from their geodesic values. We thus obtain the O(fi) correction to the specific 
energy and angular momentum parameters (in the Lorenz gauge), as well as the 0(/i) shift in the 
orbital frequency (which is gauge invariant). 
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I. INTRODUCTION AND SUMMARY 



The problem of calculating the back-reaction force, or self force (SF), experienced by a point particle as it moves 
in curved spacetime is now understood well enough to allow actual computations of this effect in systems comprising 
of a small object orbiting a large black hole. The fundamental formulation of the problem and its solution was 
set in works by DeWitt and Brehme [l[ (electromagnetic SF), Mino, Sasaki and Tanaka Q and Quinn and Wald 
3 (gravitational SF), and Quinn [i[ (scalar field SF). An alternative formulation was introduced by Detweiler and 
i Whiting [5(, also clarifying the relation between the SF picture ("forced motion on a background geometry") and 
the standard description based on the principle of equivalence ("geodesic motion in a perturbed geometry"). A 
, number of authors later devised a practical calculation method for the SF in black hole spacetimes — the "mode sum 
scheme" — which is based on multipole decomposition of the retarded field, and relies on standard methods of black 
(j [ hole perturbation theory [|| 0, H, @, G3 • This method has since been implemented by various authors on a case-by 
O" 1 case basis, so far mostly for calculations of the scalar field SF. Work so far included the cases of a static particle 
in Schwarzschild [TTf| or along the rotation axis of a Kerr black hole [l^ ]; radial plunge trajectories [l3| and circular 
bJ[). orbits around a Schwarzschild Black hole [IJ, [HI, [l6|, [I?); and ongoing work on eccentric orbits in Schwarzschild 
111 EH- The gravitational SF has been calculated so far only for radial trajectories in Schwarzschild [2(| and for 



static (non-geodesic) particles in Schwarzschild [2l| . The case of an orbiting particle has been tackled only under the 



post-Newtonian (PN) approximation [22|. A comprehensive review of the subject, including a self-contained account 
of SF fundamentals, is provided by Poisson [23|. For a snapshot of the current activity in the field, the reader may 
refer to H. 

One of the main motivations for the work on self forces draws from the need to devise accurate theoretical waveforms 
for the gravitational radiation from extreme mass-ratio inspirals (EMRIs) — of the prime targets for LISA, the planned 
space-based gravitational wave detector (25|. This requires solving the SF problem in the gravitational case, for generic 
inspiral orbits around a Kerr black hole. The main challenge in extending the analysis from the scalar-field toy model 
to the gravitational case has to do with the gauge freedom in the latter case. The problem can be summarized as 
follows. The gravitational perturbation in the vicinity of the point particle is best described using the Lorenz gauge 
(see Appendix |A"| . which preserves the local isotropic nature of the point singularity. On the other hand, the field 
equations that govern the global evolution of the metric perturbation are more tractable in gauges which comply well 
with the global symmetry of the black hole background — best known examples of which are the "radiation" gauges 
[2r| or the Regge- Wheeler gauge [13] ■ Now, in calculating the local SF we need, essentially, to subtract a suitable 
local, divergent piece of the perturbation from the full (retarded) perturbation field. In doing so, both fields (local 
and global) must be given in the same gauge; the "gauge problem" arises since the two fields are normally calculated 
in different gauges. Indeed, the only fully- worked-out example of the gravitational SF so far is the case of radial orbits 
in Schwarzschild [20| . where the gauge problem is avoided simply because, in this particular setup, the singular piece 
of the Regge- Wheeler perturbation happens to coincide with that of the Lorenz-gauge perturbation. 

One approach to the problem has been to try and calculate the local div erg e nt p iece in one of the "global" gauges — 
specifically the Regge Wheeler gauge in the Schwarzschild spacetime [111 [23, l30j |. This has been implemented so far 
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only within a post-Newtonian approximation |22j . In the current work we take a complementary approach: We solve 
the perturbation equations, and obtain the "global" retarded field, directly in the Lorenz gauge. The calculation is 
therefore done entirely within the Lorenz gauge, the "subtraction" procedure necessary for constructing the SF is 
implemented in a straightforward way, and the gauge problem is avoided altogether. Other advantages of working in 
the Lorenz gauge include the fact that the field equations then take a fully hyperbolic form (making them especially 
suitable for time-domain integration); and the fact that the Lorenz-gauge metric perturbation is better behaved near 
the particle compared with the perturbation in other gauges [3l| (which, again, makes it more suitable for numerical 
implementation) . The better reg ularity of the Lorenz-gauge perturbation is manifested in the behavior of individual 
multipole modes of the field[47|: It is well known, for example, that the multipole modes of the Regge- Wheeler 
perturbation from a point particle in Schwarzschild generally show a discontinuity across the particle. In contrast, 
the modes of the Lorenz-gauge perturbation are always continuous at the particle. 

Our "all-Lorenz-gauge" approach is made possible (at least is the Schwarzschild case) following a recent work by 
Barack and Lousto [32| (hereafter BL), which provided a practical formulation of the Lorenz-gauge perturbation 
equations in the Schwarzschild geometry. Our calculation is based on the BL formulation, and our numerical code 
incorporates the code developed in BL (with a few improvements) . In the current work we focus on circular geodesic 
orbits, for simplicity. However, since our treatment is based on a time-domain evolution, our code could be amended 
to deal with generic orbits (in Schwarzschild) in an almost straightforward manner. We shall discuss this extension 
of the analysis in the concluding section, and also comment there on the important extension to the Kerr case. 

Our calculation of the SF proceeds as follows. We first write down the 10 (coupled) evolution equations for the 
10 tensorial-harmonic components of the Lorenz-gauge metric perturbation, in the form given in BL (with a slight 
modification) . The energy- momentum of the orbiting particle is represented by a suitable delta-function distribution, 
whose tensor-harmonic components serve as sources for the evolution equations. For given orbital radius and given 
multipole numbers I and m we solve the equations numerically through time-domain evolution in 1+1 dimensions 
(time+radius), using a 2nd-order-convergent finite-difference scheme on a staggered grid based on characteristic 
coordinates. We integrate long enough to allow any spurious initial radiation to dissipate efficiently (this takes ~ 3 
orbital periods for strong-field orbits). We then record the values of the metric perturbation and it temporal and 
radial derivatives at the location of the particle. [Recall that individual multipole modes of the perturbation are 
continuous at the particle. Their first derivatives have a finite jump discontinuity across the particle (in the radial 
direction) and so we record both values of the derivatives.] We repeat this calculation for all multipole modes with 
2 < I < £maxi where Z max is determined experimentally so that our standard of accuracy (of < 10 -3 fractional error in 
the final SF) is met. In practice we found it sufficient to take Z max = 15 for the radial component and Z max = 5-9 for 
the t component (depending on the orbital radius). The modes I — 0, 1 are calculated separately, using the method 
of Detweiler and Poisson [331 ] . The values of the metric perturbation and its derivatives at the particle are then used 
as input for the mode-sum scheme. Within this scheme, each of the modes is "regularized" using functions known 
analytically [1, SOU) and the sum over modes yields the desired SF. 

Our main results are summarized in Tables ITVl and fVl (along with Fig. [8]). The tables display the values of both 
radial and temporal components of the SF as a function of the orbital radius. The temporal component of the SF 
is simply related, in our stationary circular-orbit setup, to the rate of change of the orbital energy parameter, and 
hence to the flux of energy carried in gravitational waves to null infinity and down the event horizon. Our results 
demonstrate this energy balance, which provides a reassuring validity test for our code. The radial component of the 
SF (which is itself gauge dependent) describes the conservative back-reaction effect on the orbital parameters. Based 
on our results we calculate the conservative shift in the energy and angular momentum of the circular geodesic, as 
well as the shift in the orbital frequency — the latter being gauge invariant. Our results for the shifts are plotted in 

Figs. [U and nni 

To check the validity and robustness of our code, and assess the accuracy of our results, we performed the following 
tests, (i) Numerical convergence: For each of the modes calculated, we repeated our computation with a handful of 
different numerical resolutions, checking that the answer converges quadratically to a limiting value with decreasing 
step size. To determine the limiting value we used a Richardson extrapolation (over step size), and recorded the 
estimated error from this extrapolation, (ii) Effect of spurious initial waves: We compared the values of the SF at two 
different (late) evolution times (recall that in our stationary setup the physical SF does not depend on time) , in order 
to asses the effect of residual spurious waves, (iii) Large I behavior: The behavior of the SF modes at large multipole 
numbers I is known analytically with high precision [see Eqs. (|20p - (|22|) below]. We verified that our numerically 
calculated modes have the right behavior at large I, through all three known leading terms in the l/l expansion. This 
agreement is necessary, in fact, for a successful implementation of the mode-sum scheme, (iv) Error from large-l tail: 
The mode-sum scheme involves summation over all modes I. In practice we computed all modes up to I = / max , and 
used an extrapolation to estimate the contribution from the remaining I > / max tail. We assessed and recorded the 
error from this extrapolation, (v) Comparison of one-sided forces: The mode sum scheme can be implemented in 
two essentially independent ways, by using either the "external" or the "internal" values of the SF modes (i.e., values 
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calculated by taking one-sided derivatives of the metric perturbation from outside or inside the orbit, respectively). 
Of course, the final value of the SF should not depend on our choice. We used our code to work out both values, 
and checked that they are the same to a very good accuracy. We recorded the (tiny) difference between the two 
values. Our total computation error was taken to be the combined error from the extrapolation over step size [item (i) 
above] , the deviation from stationarity [item (ii) above] , the extrapolation over I [item (iv)] , and the small discrepancy 
between external and internal values. We made sure that this combined error is kept under 0.1%. (vi) Comparison 
with energy flux: We checked that the computed temporal component of the SF balances the flux of energy to null 
infinity and down the horizon. We found an excellent agreement. 

The structure of this paper is as follows. In Sec.|TT]we review the formalism for constructing the metric perturbation 
in the Lorenz gauge and for calculating the SF via the mode sum scheme (focusing on the case of circular orbits in 
Schwarzschild). We also discuss the effect of the SF on the geodesic parameters (energy, angular momentum, angular 
velocity), and how these depend on the choice of gauge. Sec. Mil describes our numerical method, including a detailed 
description of the finite-difference scheme. In Sec. IIVI we present a few validation tests for our code, and explain how 
we estimated the computation error. Sec. |V] gives the results: We tabulate both temporal and radial components of 
the SF as a function of the orbital radius, and calculate the shift in the orbital parameters due to the conservative 
piece of the force. Finally, in Sec. I VII we discuss the extension of this work to more general orbits in Schwarzschild, 
and to orbits in Kerr. 

Throughout this work we use standard geometrized units (with c = G = 1) and metric signature ( — | — h+). The 
Riemann tensor is defined as in Ref. [34| . 



II. REVIEW OF THEORY: SELF FORCE IN LORENZ GAUGE 



A. Orbital setup and equation of motion 



Consider a pointlike particle of mass fi, in a circular orbit around a Schwarzschild black hole with mass M 3> /i. 
Let the worldline of the particle be represented by x a — x^{t), with tangent four velocity u a — dx^/dr. At the 
limit jti — * (i.e., neglecting SF effects) the particle traces a geodesic xp* = Xq(t), with associated four velocity 
Uq = dx^/dr. Without limiting the generality, we adopt a Schwarzschild coordinate system t,r,6,tp in which the 
orbit is confined to the equatorial plane. Then 

x o ( T ) = [*o(t), r = const, 6 = vr/2, <£>o(t)] . (1) 

This circular geodesic can be parametrized by the radius r"o, or, alternatively, by the angular velocity (with respect 
to time t) 

il Q = dipo/dt = y/M/r$, (2) 

by the specific energy parameter, 

£ ee -u ot = (1 - 2M/r )(l - 3M/r )- 1/2 , (3) 

or by the specific angular momentum parameter, 

C ee u Qv = (Mr ) 1/2 (1 - 3M/r )- 1/2 . (4) 

The subscripts '0' here indicate that the above are values associated with the geodesic Xg (below we will consider the 
correction to these values due the SF effect). As always in our perturbative treatment, tensorial indices are "raised" 
and "lowered" using the background metric. 

Now assume that /x is finite (but still very small compared to M) . The equation of motion can be written as 

D 2 x" Du a 

where the covariant derivatives are taken with respect to the background (Schwarzschild) geometry, and _F"[~ 0(^ 2 )] 
is the gravitational SF. Clearly, the symmetry of the problem implies F e — 0. Also, assuming the four-velocity is 
kept normalized along the worldline, i.e., u a u a = — 1 , we have D(u a u a ) / Dt = 0, leading to u a F a — 0, and the four 
components of the SF are not independent. In the circular orbit case we have the relation UtF + u^F^ — 0, which 
we may write, through leading order in fi, as 

F v - {Bo/C^F*. (6) 
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Hence, we need only calculate two of the components of the SF: the r component and (say) the t component. For 
simplicity we shall refer to these as the "radial" and "temporal" components. Our goal would be to calculate both 
components, as a function of the orbital radius ro- Note that it is sufficient, for the sake of obtaining the leading-order 
[0(/i 2 )] SF, to assume that the motion is momentarily geodesic. 

The SF affects the motion of the particle in two ways: Firstly, —u t and u v are no longer conserved over time, so 
we can speak of the "rate of change" of the energy and angular momentum of the orbit. Secondly, at each given 
time, the values of —u t and u v are shifted with respect to their corresponding geodesic values — Mot and uq v . In the 
case of a circular orbit, the first, "dissipative" effect is due entirely to F t (and F v ), while the second, "conservative" 



effect is due entirely to F r 
immediately give 



To see this, start by defining £ = —ut and C 



The t and ip components of Eq. 



If 



dC 



= -if 1 Ft and — = fi^E 



dr 



(7) 



respectively, describing the dissipative effect of the SF. The change of energy and angular momentum is precisely 
balanced by the flux of energy and angular momentum carried away by gravitational waves. For the conservative 
effect, consider the r component of Eq. along with the normalization condition u a u a = — 1. Solving these two 
equations simultaneously for u* and if (Recalling u r = du r /dr = 0), one obtains (it 4 ) 2 = ro (1 — roF r /n) /(ro — 3M) 
and (u v ) 2 = (M/rl - F r /fj) /(r - 3M), or, through 0(fi), 



£ = £o 



F r 



C 



( rl 



F' 



(8) 



Given the radial component of the SF, the last two equations give the "conservative" shift in the energy and angular 
momentum parameters. It is also useful to look at the shift in the orbital frequency fl = dip p /dt = u^/u 1 . Based on 
the above expressions for u* and u v we obtain, through 0(/x), 



i 



r (r - 3M) 
2M/i 



F r 



(9) 



which describes the shift in the "frequency at infinity" due to the conservative piece of the SF. 



B. Gauge dependence 



It is important to understand how the above quantities depend on the choice of gauge. Let ft. a| a[~ O(fi)] be the 
metric perturbation due to the particle, given in a specific gauge. Consider an infinitesimal gauge transformation 



This will change the metric perturbation by an amount 

d^hap = h a p — h a p 
It will also induce a change in the SF, given by [3l| 



-(£a;/3 + £/3;a)- 



5 ( F a =F' a -F a = l x 



[9a\(xo] + U 0a U x 



2tA 



Dt 2 



(10) 



(11) 



(12) 



Here g a \{xo) and R ai i,\v(xo) are the background (Schwarzschild) metric and Riemann tensors, respectively, evaluated 
at the particle. 

We now assume that the original gauge is "physically reasonable". In our case (a circular equatorial orbit in 
Schwarzschild), this would mean that the perturbation h a p in that gauge reflects the stationarity of the problem, 
and also retains the symmetry of reflection through the equatorial plane. The Lorenz gauge (see Appendix [X] for 
definition) is an example of such a gauge. Restricting our discussion to gauge transformations within the family of 
"physically reasonable" gauges, we should clearly require, in our case, d£ a /dr — 0, as well as £ e = 0. From Eq. (fTj 
we then get S^F t = 5^F$ = S^F V = 0, along with 



(13) 
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Hence, "physically reasonable" gauge transformations may affect the radial component of the SF (they do so if they 
have £ r ^ 0), but not the other components. 

One immediate consequence of the above is that the quantities d£/d,T and dC/dr in Eq. (|7|) are invariant under a 
gauge transformation (as should be expected on physical grounds). However, the quantities £ and C themselves are 
not gauge invariant. To see this, use Eq. (fT3|) in Eq. (J8]), along with £ — > £o + {d£o/dro)C , and Cq — » Co + (d£o/dro)£ r . 
[Note here that the coordinate location tq (and hence also £o and Co) is obviously not gauge invariant: Under the 
gauge transformation (|10p we have simply tq — > ro + £ r .] Denoting ^£ = £' — £ and = £ — C, and keeping only 
terms linear in £ r , one then obtains the following gauge transformation formulas for £ and C: 

s £= _2M/^ = 2 

VI - 3A//r VI - 3M/r„ 

We can construct orbital parameters that are invariant under the transformation (|10[) . One example is the frequency 
f2, given in Eq. ©. Detweiler pointed out recently [35| that the combination £ — D.C = S is also gauge invariant (for 
circular orbits). For both SI and S, it is a straightforward exercise to show that, under the gauge transformation (|10p . 

5 s O = 0, <5 S £ = 5 S {£ - QC) = 0. (15) 

C. Metric Perturbation in Lorenz gauge 

In this work we calculate the SF in the Lorenz gauge. This will require knowledge of the full (retarded) metric 
perturbation h a p in the Lorenz gauge. We briefly review here the construction of the Lorenz-gauge perturbation, 
referring the reader to BL [32] for further details. 

In the formulation by BL, the Lorenz-gauge metric perturbation in Schwarzschild is constructed through [48| 

oo / 



2r 

1=0 m 



with 



h\™ = (hW +fh (6 AY lm , 

h l ™ = r 2 (h {1) - Y l 
h l $ = r (h^Y^ + Py^ 1 



h$ = rsm0 



n r8 — 

7 Im 
ll rtp 



Urn _ 2 



rf- 1 (h&Y#? + h i9) Y^ , 
rf- 1 sin8 (h^Y<% - h^ 9] Y^ , 



h l $ = r 2 sind (h^Y^ -h^Ytt 



h%l = r 2 sin 2 i 



Umylm _ fc(7)yfan _ h (W) Y ±™\ . (17) 



Here / = 1 — 2M/r, and ly™, Yy%, Y^ , and Y^,™ are angular functions constructed from the standard spherical 
harmonics Y lm (9, ip) through 

YiK = j^jj Yf (forZ>0), 

r v™ = J(iTT) sin_1 6Y « (for Z > 0) ' 



\rlm 

r T1 



[sin0 (sin- 1 0y<, m ) )e - sin- 2 y'™] (for I > 1), 
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y ^ = jtw ^^y.e (for 1 > i} - (18) 

The functions h^ lm [i — 1, 10; the indices I, m were omitted in Eq. (| 1 7[) for brevity] depend on r and t only, 
and form our basic set of perturbation fields. These time-radial fields are obtained as solutions to the coupled set of 
hyperbolic [in a 2-dimcnsional (2-D) sense], scalar-like equations 

nh® lm + M^h® 1 ™ = S® lm (i= 1, . . . , 10). (19) 

Here '□' represents the 2-D scalar field operator, □ = d uv 4- V(r), where v and u are the standard Eddington- 
Finkelstein null coordinates, and the potential is V(r) = \f [2M/r 3 + 1(1 + l)/r 2 ] . The "coupling" terms M%h^ lm 

involve first derivatives of the h^ lm, s at most (no second derivatives), and gM'™ are source terms. Both Ai )^.h^ lm 

and gW™ are given explicitly in Appendix [Al for our circular orbit case. In addition to the evolution equations (119p . 
the functions ftM im also satisfy four elliptic equations, which stem from the gauge conditions. These "constraint" 
equations are also given in Appendix El The set (fT9| incorporates "divergence dissipatin g" t erms, which guarantee 
that violations of the Lorenz gauge conditions are efficiently damped during the evolution [32J . 

BL demonstrated in [32| how Eqs. (fTi?|) can be evolved numerically (for I > 2), in the time domain, with a delta- 
function source term representing an orbiting point particle. They also demonstrated that the solutions preserve the 
Lorenz gauge condition throughout the late-time evolution (initial violations of the gauge conditions are suppressed 
over a time- scale ~ M) . 

The perturbation modes I = and I = 1 require a separate treatment. For I = 0, the set (ITO)) reduces to 4 equations 
(for ft/ 1,2,3,6 )), which describe the spherically-symmetric, monopole mass perturbation. Sec. III-D of BL gives the 
solution for the Lorenz-gauge monopole perturbation in analytic form (based on analysis by Detweiler and Poisson 
|33l]). For I = l,m= ±1, the set ((THJ) reduces to 6 equations (for h^ 1 ^ 6 - ) ). These describe the rotational dipole piece 
of the perturbation, which is (in Newtonian terms) due to the motion of the black hole about the center of mass of 
the black hole-particle system. The Lorenz-gauge solution for this mode was obtained by Detweiler and Poisson [33[ , 
using a procedure that reduces the problem to the solution of 3 coupled ordinary differential equations (see also Ori 
[36| for a fully analytic weak-field treatment). Finally, for I = l,m = — the axisymmetric dipole perturbation — the 
set (|19p reduces to a single equation (for h^), describing the perturbation in the angular momentum due to the 
particle. Analytic solution for this mode was obtained long ago be Zerilli [37|], and identified in [33] as a Lorenz gauge 
solution. BL later obtained analytic Lorenz-gauge solutions for all axisymmetric (m = 0) modes with odd I. These 
are given in Sec. III-C of BL. 

In what follows we will prescribe the construction of the gravitational SF directly in terms of the functions h^ lm . 



D. Self force via the mode-sum method 



In the mode sum scheme, the gravitational SF is calculated through 0, d, H E3] 

00 



1=0 



(20) 



where L = 1 + 1/2, and Fgi, are the modes of the "full" force, obtained from the modes of the full (retarded) 
metric perturbation in a manner described below. The subscript ± refers to the two possible values of at xq, 



resulting from taking one-sided radial derivatives of the metric perturbation from cither r — > 7q or r — > r$ (recall that 
the modes of the Lorenz-gauge metric perturbation are continuous, but their gradients generally have a finite jump 
discontinuity across the particle). A± and B a are the "regularization parameters", which are known analytically. For 
circular equatorial geodesies in Schwarzschild we have A± = B a = for a = t,9, <p, and 



A'± = =F*r ( 1 - — 



(21) 



D' 



n 2 F 2 



^+r§)3/2 



E(w) - 2K(w) 



(22) 



where K(w) = J (1 — wsm 2 x) x l 2 dx and E(w) = J (1 — w sin 2 xY^dx are complete elliptic integrals of the 
first and second kind, respectively, and w = (tq/M — 2) _1 . The final SF F a can be calculated using either "external" 
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(+) or "internal" (— ) values; The difference [i ? ftji(aro)] ± — A±L is guaranteed to be direction-independent. One can 
also use the average value, Fjgh = {[F^ n ] + + [F^ u ]-} /2, in terms of which the mode -sum formula takes the more 
compact form 

oo 

F a = £ [fgjxCro) - B a ] . (23) 
;=o 

Recall here we are interested in the r and t components of the SF. For the r component, the contribution from the 
individual full modes is oc L at large i, and the sum over [-Ff„u(xo)] ± diverges. However, the contribution from the 
"regularized" modes [-Pfuu^o)] ± — A r ±L — B r falls off as oc L~ 2 at large L, and their sum converges (and gives the 
correct physical SF). The regularized modes admit the large-Z expansion 

F;4 = [F[l n (x )] ± -A r ± L-B r = ^ + ^ + ---, (24) 

where D2, -D4, . . . are coefficients that may depend on the orbital parameters, but not on L. As for the t component: 
In the special case of a circular orbit the temporal full- force modes require no regularization (recall A t — B l =0), 
and their sum converges. In fact, in this case the mode sum can be shown to converge exponentially at large I (this 
will be demonstrated experimentally in Sec. IIVI below). For later convenience we write 

and the mode-sum formula becomes, for cither the r or the t component, 

00 

F a = Y, F rt- (26) 



E. Construction of the full-force modes 

We now describe the construction of the full modes [i^iO^o)] + out of the Lorenz-gauge perturbation fields h^ lm . 
First, following 0, define the "full force field" as a tensor field at arbitrary spacetime point x, for a given (fixed) 
worldlinc point xo (where the SF is to be calculated): 

F£ u (z;z )=m^ t5 W- (27) 
Here h a p = h a p — \g a pg^ v h^ v is the "trace-reversed" Lorenz-gauge metric perturbation at x, and 

k a ^\x;x ) = g a5 u p u'/2 - g ap u*<u s - + u a g^u s /4, + g aS g M /A, (28) 

where g aS is the background metric at x, and u a are the values of the contravariant components of the four- velocity 
at £0 (treated as fixed coefficients). Obviously, the full force F f " n diverges for x — > xq (like ^distance -2 ). 

Next, expand the metric perturbation in tensor harmonics as in Eq. (|16|) . and substitute in Eq. I|27j) . Taking the 
limits r — > ro and t — » to (but maintaining the 9, tp dependence), the full force takes the form 

2 00 1 

[f£u(0,¥>;r o ,to)] ± = E {/o± m ^ m + fit sin 2 6Y lm + f& m co S sin 9Y$> 

1=0 m=-l 

+ f^ n sin 2 6 Y% 1 + ff ± (cos 9Y lm - sin 8Y l g m ) 

+ f^ n sin 6 Y 1 ^ + f^ n sin 3 6 Y l e m + f^ m cos 6 sin 2 6 Y$ } , (29) 

where Y lm (6,tp) are the spherical harmonics, and the coefficients f%± are constructed from the perturbation fields 
Jii l ) lm anc j their first r and t derivatives, all evaluated at xq- The labels +/— correspond to taking external/internal 
r derivatives, respectively. The explicit expressions for the /"j_ m 's are quite lengthy, and we give them separately, in 
Appendix IB"1 

The individual I modes in Eq. (|2"9"|) are not quite yet the full force modes needed in the mode-sum formula (|20[) : 
A little complication arises because the mode-sum formula requires the decomposition of the full force in scalar 
harmonics. That is, to obtain the modes [-ffun( x o)] ± f° r use m Eq. (|20p we are required to ignore the vectorial nature 
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of the full force, and expand each of its components in scalar harmonics. To obtain this, we need only re-expand all 
angular functions in Eq. (|29[) in spherical harmonics, and rearrange the terms in the sum. The necessary expansion 
formulas are given in Appendix O We find that each of the tensor-harmonic I modes in Eq. ([29]) couples to a finite 
number of scalar harmonics (there is no coupling between different m modes). For the r component, each tensor- 
harmonic I generally couples to 5 scalar harmonic modes V with I — 2 < V < I + 2; for the t component, it generally 
couples to the 7 modes / — 3 < I' < I + 3. It should be commented that this finite coupling is characteristic of all 
(equatorial) orbits in Schwarzschild, not necessarily circular (although in general both t and r components would 
involve coupling to 7 scalar harmonics). We also comment that, in our circular orbit case, it is not quite necessary to 
expand the t component in scalar harmonics, since this component requires no regularization and so the mode-sum 
can be evaluated directly using the tensor harmonic expansion. In this work we nevertheless choose to expand the 
t component too in spherical harmonics, for two reasons: Firstly, this would allow us to test our general treatment 
of the coupling between modes, since the computed value of the t component of the SF could readily be verified 
by comparing with the flux of radiated energy. Secondly, our code could then be more easily adapted to deal with 
eccentric orbits, where scalar-harmonic decomposition of both r and t components is necessary. 

Hence, we now re-expand Eq. (|29|) in spherical harmonics (using the relations given in Appendix [C]), then rearrange 
the terms in the sum by collecting all terms with the same scalar-harmonic multipole number I, and, finally, set 9 — > 9q 
and tp —> ipo. The resulting "scalar harmonic" / modes of the full force take the final form 

2 * 

(„. W ^ \ y \rlm(n ,„ \ { T-ai-3,m . ^-al—2,m . ~~al — \,m . -raZm , -j-ctZ+l.rn . ^-al+2.m . ral+3,m\ 

L^fullC^OjJ ± - -J Y WW X ^ 3 ) +J (-2) + +J (+2) L' 

T m=-l * 

(30) 

where 



■>-{-S) 

aim 



(-1) 



(0) 



<r-alm 

<r-alm 
A+3) 



Clm ralm . 1 1 m ralm 
;+3)/6± " t "^(+3)/7± ) 

„ Im fOtXm I aim .calm t _,/m r a I m 
tt( + 2)/l± +P(+2)/2± +7(+2)/3± J 

Im ralm , rim palm , Am ralm , Am ralm 
e (+l)J4± -r °(+l) /5± "T S(+l)/6± "T S(+l)/7± 

r a I m , „ aim, ralm . olm £alm , „,lm ralm 

Jo± + a (o) fi± +P(o)h± +7( )/3± » 

Im ralm . clm ralm ■ Am ralm ■ Am r 
6 (-l)/4± "T J5± ' S(_l)/6± ' 

„ Im ralm , olm ralm , „Jm ralm 
a (-2)Jl± +P(-2)/2± +7(-2)/3± ! 

Am falrn 
>(-3)/7± ' 



<\trn 
7± ! 



S.(-3)/6± 



(31) 



with the various coefficients a, /3, 7, (5, e, £, and £ given in Appendix [Cl Note /g ± = /f ± = 0, so that Eq. (f30f simplifies 
considerably for the r component. The spherical harmonic Y lm (8o,(po) is given explicitly, for 8q = tt/2, by 



y im (7T/2^ ) 




(2Z+l)(i+-m-l)!!(i-m-l)!! 
4ir({+m)!l(i-m)!l 



1/2 



/ — m even 
/ — m odd. 



(32) 



Hence, for given I, only m modes with even / — to contribute to the sum in Eq. (|30[) . A further simplification arises 
since the individual to modes in the sum in Eq. (|30[) arc invariant under to — > —to, which allows us to fold the to < 
part of the sum over to to > 0. In practice, therefore, one is required to compute only 2/2+1 TO-modes for each even 



2-mode, and (I + l)/2 TO-modes for each odd 2-mode. Finally, note a^ 2 ^ = /3[+ 2 ) = 7(+2) = e (+i) 
= Cf+31 = = for I < 0, such that no functions f% lm with I < occur in Eq. J3DJ). 



= (5 



(+1) 



(+1) 



F. Summary of prescription for constructing the Lorenz-gauge SF 

Start by calculating the Lorenz-gauge perturbation fields h^ Lm (10 of which for each I > 2, to), by solving (numer- 
ically) the field equations (p~9|) . Obtain the modes I = 0, 1 of h^ lm through the procedure described at the end of 
Sec. Ill CI Construct the functions and using the formulas in Appendix [Bl and then construct the quantities 
T through Eqs. (f3"Tj) . Use these in Eq. (f3TJ)) to obtain the scalar-harmonic 2-modes of the full force, [-PSu(^o)] ± - 
Incorporate the full force modes in the mode-sum formula (|20p [or (|23[) ] to obtain the SF. 
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III. NUMERICAL IMPLEMENTATION 



In this section we summarize the numerical method used to calculate the SF through the mode-sum formula (|20[) . 
The evolution of the Lorenz-gauge field equations (fill)) was described in BL, but we will briefly review this method 
here — mainly in order to supplement a few details of the finite difference scheme left out in BL. We then describe 
the numerical construction of the regularized SF modes, and the technique used to evaluated the infinite sum over I. 
There are three main sources of numerical error in our calculation: (i) error from the finite-grid discretization (which, 
in the procedure described below, comes from the error in a Richardson- type extrapolation to zero grid size); (ii) error 
from estimation of the large-/ tail of the mode-sum series; and (iii) error from residual spurious waves resulting from 
the imperfection of initial data. We explain how all these errors are monitored and controlled in our calculation. 



To solve Eqs. ([19]) for the various modes I > 2 we use characteristic time-domain evolution on a fixed 2-D staggered 
double-null grid based on v, u coordinates. The numerical domain is depicted in Fig. 2 of BL. The evolution starts 
with characteristic initial data on two initial "rays" v = va and u — uo, taken such that the vortex vq, uo corresponds 
to r = ro with initial time to = 0. The circular orbit worldlinc then traces a straight vertical line through the grid, 
connecting the "lower" and "upper" vertices (see Fig. 2 of BL). In this setup the worldline cuts through grid points. 
This does not cause problem, since the Lorenz-gauge perturbation modes are continuous at the worldline. For initial 
data we simply take h^ lm = along v — vq and u — uq, for all i. This sparks a burst of spurious radiation at the 
initial vortex, which, however, dies off efficiently over time and leaves very little trace after 1-2 orbital periods of 
evolution (we demonstrate this below). 

Our finite-difference scheme (particularly the handling of the delta function source term) is based on the method first 
introduced by Lousto and Price [38| and later implemented by a number of authors [HI, [2(| HjJ [4(| • In this method, 
the finite-difference equation is obtained by approximating the (2-D) integral of both sides of the field equation over 
a grid cell, at a suitable accuracy. This automatically deals with the delta-function singularity on the right-hand side 
of the equation. For the following discussion, consider Fig. [T] Suppose that we have already solved for all h^' lm, s at 
the grid points numbered 2, 3, and 4. Denote the values calculated at these points by h% \ , and , respectively 
(we omit here the indices l,m for brevity), and let the sides of the grid cell be Av — Au — h. We are interested in 

obtaining h^, the value of the h^'s at point 1. 



FIG. 1: A numerical grid cell, of dimensions h x h (see description in the text). Our 2-D grid is based on characteristic 
(Eddington-Finkelstein) coordinates v and u. These are related to the Schwarzschild coordinates through v = t + r* and 
u = t- r», where r» = r + 2M\n[r/(2M) - 1]. 

Consider first the principal part of the field equations (fl9|) . For any of the ten i's, it reads h\uv- This term is 
integrated exactly over the grid cell, to give 



A. Metric perturbation: finite-difference scheme for I > 2 



1 





(33) 



Since the h^>'s are continuous at the worldline, the above integral holds even for grid cells crossed by the particle. 
The remaining part of the left-hand side of the field equations p^|) includes three types of terms, of the form V\(r)h^ , 
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^(rj/i'r , and Vs{r)h/^ , where the \^(r)'s are known radial functions. As for terms of the first two types, we can 
approximate their integrals over the grid cell as 

/ v^m ** . i*nw<*s« + *> + { gj$ <«> 

/ VMHf ** . ^rJWj? - 5<") + { g<$ ™ le)i (35) 

where r c is the value of r at point C in the middle of the cell (see Fig. QJ. The case indicated as "particle" refers 
to grid cells crossed by the worldline. "No particle" refers to grid cells not crossed by the worldline. The difference 
in the error terms arises because the ftW's generally have discontinues r derivatives across the worldline. We can 
accommodate a local error term of 0(h 3 ) along the worldline, as the worldline is crossed only once in each evolution 
time step, and so such an error accumulates over time to give a global error of only 0(h ). 




Case A Case B 



FIG. 2: Finite-difference scheme for terms in the field equations involving single v derivatives (diagrams to go with the 
description in the text). The dashed line represents the worldline, with two cases shown. 

Terms of the form V 3 (r)h!if require a more careful treatment. For these we will need information from outside the 
grid cell of Fig. Q] Consider Fig. [2l showing an extended area around the central point C, now including also points 

- U\ - u\ - n\ 

5-9. We assume all functions h 2 h g have been calculated before, and we need to obtain h\ . The figure shows two 
special cases: In "case A" the worldline crosses the point C; in "case B" it crosses "just to the right" of point C, 
through points 3 and 6. The integral over the term V^(r)h\v can be approximated, in the various cases, through 

f 2hf - hf - fhf - \hf + 0(h 3 ) (Case A), 

V 3 (r)hf dudv = -hV 3 (r c ) x } hf - hf + 2,'hf - 2hf - 2hf + hf + 0(h 3 ) (Case B), (36) 
cc " " [ 3hf - 3hf - hf + hf + 0{h 4 ) (no particle), 

where "no particle" refers to the most common case, in which the particle passes through neither point C nor point 
3. There are, of course, alternative schemes which approximate this integral at the same order of accuracy. We have 
chosen the particular scheme in (|36p as we found it experimentally stable. 

Finally, we need to integrate the source term S^ lm on the right-hand side of the field equations [the explicit form 
of the source term is given in Eq. (|A11[) of Appendix [A] . Thanks to the delta function in S^' we can work out the 
integral over the grid cell exactly. We find 

S {l)lm dudv = 8ir£ a {t} r 1 (r c )xhxsinc(mn a h/2) 

ill 

y p -imn t a y f Y lm *(n/2, 0), for i = 1-7, 

Xe X \ y;(7r/2,0), fori = 8-10, [dt) 

where t c is the value of t at point C, an asterix denotes complex conjugation, the coefficients a' 1 ' are those given in 
Eqs. (|A12|) . and sinc(x) = (smx)/x for any i^O, with sinc(0) = 1. 
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Integrating the field equations (fT9|) using above Eqs. (|33)) - (f37|) . we can solve for the ft^'s at point 1 given the 
values calculated in previous steps of the evolution. In this scheme we have to keep two w=const data vectors at 
each (advanced) time step. The local finite-differentiating error at each grid point scales as ~ ft 4 , except for points 
belonging to Cases A or B (Fig. [5]), for which the local error scales as ~ ft 3 . Since the total number of steps scales as 
~ ft~ 2 , and the number of steps belonging to Cases A or B scales as ft- -1 , we expect the error accumulated over the 
entire evolution to scale as ~ ft 2 . 



B. Metric perturbation: Monopole and dipole modes 

The monopole and dipole modes are dealt with separately. For I = we use the analytic solution for the ftW's 
from Sec. III-D of BL. For the mode I = 1, m = we use the analytic solution from Sec. III-C therein. For the mode 
I = m = 1 we follow the method of Detweiler and Poisson [33| , which involves solving (numerically) a coupled set 
of 3 ordinary differential equations, with boundary conditions at infinity and along the horizon, and with matching 
conditions across the particle. 

One may attempt to compute the modes I = 0, 1 too using the evolution equations (|19[) (which reduce to four 
equations for I = 0, six equations for I = m = 1, and one equation for I = 1, m = 0). In practice, however, the 
system (fl~9|) . in its present form, does not seem to evolve stably for these modes: For I < 2, some of the potential 
functions in these equations turn negative for some r values outside the black hole, apparently rendering the evolution 
unstable. This is not a serious problem in our present analysis, as we simply derive these two modes using other 
methods. The problem will have to be address when extending the analysis to non-circular orbits, for which analytic 
(or semi-analytic) solutions are not yet at hand. One may then either attempt to derive analytic (or semi-analytic) 
I = 0, 1 solutions for eccentric orbits; or, alternatively, attempt to find a re-formulation of the evolution equations 
suitable for I = 0, 1. 



C. Taking derivatives of the metric perturbation 



- u\ - u\ 

The construction formulas for the r and t components of the SF require the derivatives hy and ft t 



at the particle. For the t derivatives we can simply make the substitution <9 t 
case the fields ft^ are stationary and depend upon t only through [ e lm voW]* c 
numerically, using the finite-difference formula 



iQnt 



3ft, 



(0 



4ft 



±i 



:(0 

l ±2 



2ft/(r ) 



+ 0(ft 2 ), 



both evaluated 
imtto, since in the circular orbit 
The r derivatives are taken 



(38) 



where the various quantities refer to the diagram in Fig. [3] ftg is the value calculated on the worldline, at a given 
time, and h±i and ft±2 are the values at points ±1 and ±2, respectively. Recall /(ro) = 1 — 2M/r$. The subscripts 
+/— refer to one-sided derivatives taken from Tq or r$. This scheme allows for discontinues r derivatives, which we 
expect some of the 's to have. The scheme gives the local derivative with an error that scales as ~ ft 2 — the same 
as the accumulated error in the TS l \ Hence, we expect the total finite-differentiation error in the SF to scale as ~ ft 2 . 




t=const. 



FIG. 3: Diagram to explain how r derivatives are taken at the worldline (see description in the text). The dashed line represents 
the particle's trajectory on the numerical grid. The SF is calculated at the point labeled '0'. 



Once we have at hand the various ftW's and their derivatives at the particle, we can construct the various scalar- 
harmonic modes of the full force, [-Ff^ n (a;o)]± and [F[^ u (xq)]±, using Eq. (f3"01) . [Recall that to obtain a single scalar- 
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harmonic mode of the full force we need to calculate all tensor-harmonic modes of the perturbation with multipole 
numbers between I — 2 and 1 + 2 (I — 3 to I + 3 for the t component).] We then construct the regularized modes 
through Eqs. ([24| and ([25| . We use this procedure to calculate all modes up to I = 15 for the radial component, 
and up to I = 4-9 (depending on ro) for the temporal component. (For the two radii ro = 6M and tq = 100M we 
obtain all modes of the radial component up to I = 25; we use the extra mode information in testing the validity of 
our numerical procedure — see below.) 

Of course, the stationarity of our problem allows us to choose any point along the orbit for calculating the force. 
We need to make sure, though, that this point is taken late enough in the evolution, where the effect of initial spurious 
waves is negligible. To monitor any residual effect from the initial waves we repeat the calculation at two different 
evolution times. We will give more details on this procedure in Sec. IIVBI below. 



D. Extrapolation to zero step size 

To obtain Fg* g with good accuracy requires to calculate [^f"n(^o)]± with an even better accuracy, which, in turn, 
requires to evolve the field equations with a sufficiently fine grid. This can become very demanding computationally, 
especially for large I. To reduce computational cost we use a Richardson extrapolation to h — > 0. The idea is to 
extrapolate the value of (using a rational function) to the limit of vanishing step size, using a sequence of values 
obtained with progressively decreasing step sizes. Specifically, we employ the Bulirsch-Stoer method (4l[, which 
utilizes the sequence 

fH = —, (39) 



ii, 



where i = 1, 2, 3, . . . and 



= {2, 4, 6, 8, 12,-"} (ni = 2n<_ 2 ). (40) 



Namely, we repeat the calculation of F£} g for all step sizes hi, h%, . . . , ft-i max , and then extrapolate the resulting series 
of values to i — > oo (h — > 0). To control the error in this procedure we introduce the estimator 



A Q 



max. 



pal \j 1 _ pal r ■ _ 1 1 
rcg L maxj reg L max 



(41) 



where F^, g [i m&x ] is the value extrapolated from the sequence of z max values of F^ obtained with hi,..., h inm . We 
repeat our calculation with increasing values of i max , until A" z [? max ] is smaller than a prescribed threshold, A^ resh . 

What value should AJ^ resh be set to? This requires some consideration. Since the contribution to the mode sum 
from individual I modes decreases with I (at large /), it makes sense to relax the threshold for large I modes. This 
is certainly true for the t component, for which the mode sum converges exponentially at large I. Accordingly, for 
the t component we set the tight threshold A*( lrcsh = 10~ 4 for each of the modes I < 3, but for I > 3 we set 



Atl _ in- 
^thrcsh — 1U 



M-l 



This is slightly larger than 10 for I = 4, but grows exponentially at large 

I. Experimentally, it yields AfL , ~ 1 for / = 5-9 (depending on tq: higher I for smaller ro). In fact, for the t 
component we use A*J lrcsh also as an indicator to tell us when it is appropriate to terminate the sum over modes: 
We sum up to the first / mode for which A*J lrcsh > 1. This guarantees an overall truncation error less than a few 
xl0~ 4 in the t component of the SF. As for the r component: Here the mode sum converges only as ~ l~ 2 , and 
determining A^ resh requires more caution. As we describe below, we estimate the contribution from the truncated 
I > 15 tail by extrapolating the numerical data from I < 15; Assigning an /-dependent threshold could make it difficult 
to control the error from such an estimation. For the r component we therefore conservatively set a fixed threshold 
of A^ rcsh = 10 -2 for each individual I mode computed. 

For both r and t components, we estimate the overall (fractional) discretization error in the SF as 



discr 





A ai [thresh] -f^og [^thresh] 




X/i=0 ^rog [*thrcsh] 





(42) 



where ithrcsh (depending on I and r ) is the smallest i max for which A Q '[z max ] < A"^ rosh . Note that we take here the 
fractional total error as the sum of fractional errors from the individual /-modes, rather than their root-mean-square 
value. This makes sense because the individual errors are not randomly distributed but rather reflect a systematic 
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extrapolation error. In practice, to reach the above thresholds, we needed ithresh values of up to 8 for the r component, 
and up to 9 for the t component (depending on / and r : larger ithresh is generally required for larger / and smaller 
to). 

For each value of r$ and / we started the above procedure with i max = 4 for the r component and i max = 3 for the 
t component; namely, we used at least 4(3) terms for the Richardson extrapolation. For some of the low-/ modes this 
already yielded A al [i max ] smaller than A^ rcsh (this is partly because the low-/ modes have large I — 0, 1 tensorial- 
harmonic components, which are available analytically). Since the total error A^ iscr is dominated by errors from these 
low-/ modes, we generally find A^ iscr values much smaller than the above set thresholds. The actual values obtained 
for Aj iscr will be stated in the Results section. 



E. Estimation of contribution from large-/ tail 

The mode-sum formula requires summation over all regularized modes from I = to I = oo. In practice, of 
course, we calculate only a finite number of modes, < I < / max . It is then necessary to estimate the contribution 
from the remaining, truncated part of the series. This is straightforward in the case of the t component, where the 
mode sum converges exponentially, and thus the contribution from the truncated tail drops exponentially with / max . 
We find experimentally that it is sufficient to take / max = 4-9 (depending on r ; larger / max is needed for smaller ro) 
for the contribution from the truncated tail to drop below our standard discretization error (~ 10~ 4 ). 

The situation is different with the r component, where the mode sum converges as cx l~ 2 , and the contribution 
from the truncated high-/ tail scales as 1/Z max . For computationally realistic values of / max (here we take / max = 15) 
this contribution cannot be neglected. To evaluate it we apply the following strategy. Let 

fr = ^+ f l>U. ( 43 ) 

where 

[ max oo 

F kl^ = T, F re B and Ff>l m = E <r ( 44 ) 



(=0 l=l„ 



The part F[ Kl is computed numerically. To evaluate Ff >l we extrapolate the last few modes in Ff^ using 
the fitting formula 



TV 



F rl ^^£2n , 45 s 
log — Z^i Tin K °l 



= 1 



(for some N > 1), where, recall, L = I + 1/2, and D r 2n are /-independent coefficients, which serve here as fitting 
parameters. In practice we have used the last 6 modes of F[ <t (i.e., 10 < I < 15) for the fitting, but have checked 
that fitting using a different number of modes does not change the result significantly (we demonstrate this below). 
Given the coefficients F) r 2n , the large-/ tail is estimated as 



Ff >lmax ~J2 D *n E ^ 2 " = E7^ nj *(2n-l,/ max + 3/2), (46) 

71=1 i = i max + l 71=1 ^ 

where ty(n,x) is the polygamma function of order n, defined as 

.(„,„_ qterfea, (47) 

in which T(x) is the standard gamma function. 

To determine how many terms it is necessary to include in the fitting formula (145[) requires some experimentation. 
The data in Tables [J demonstrate the effect of varying N: It shows the values obtained for F[ >1 (both external and 
internal values) using TV in the range 1-4. We display data for the two sample radii ro = 6M and ro = 100A/. We find 
that taking N = 3 or N = 4 gives the same value of -Fj>j as taking N = 2, with a fractional difference of merely 
< 10 at most. Since F[ >t itself contributes at most ~ 2% of the total force (see Tables IVIIVIII in Appendix [D|l . 
we conclude that it is sufficient to take N = 2. Taking only N — 1 would produce a fitting error similar in magnitude 
to Aj iscr ; So, taking N = 2 effectively eliminates the large-/ fitting as a source of error in our calculation. We find 
similar numbers for other values of ro, and so we take N = 2 in all cases. 
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[Ff >15 L 


Relative difference 


[W>is] + 


Relative difference 


N 


(fit using 10 < I < 15) 


w.r.t N = 2 


(fit using 10 < I < 15) 


w.r.t N = 2 


r = 6M 


1 


-5.117158 [2 x 1(T 3 ] 


-1.4 x 10" 2 


-5.117176 [2 x 10" 3 ] 


-1.4 x 10" 2 


2 


-5.046144 [3 x 10 -5 ] 





-5.046189 [3 x 10" 5 ] 





3 


-5.046984 [6 x 10~ 5 ] 


-1.7 x 10 -4 


-5.046812 [2 x 10" 4 ] 


-1.2 x 10" 4 


4 


-5.046316 [4 x 10 -4 ] 


-3.4 x 10" 5 


-5.046133 [9 x 10 -4 ] 


-1.1 x 10~ 5 


r = 100M 


1 


-5.904165 [1 x 10" 3 ] 


-7.5 x 10~ 3 


-5.904159 [1 x 10" 3 ] 


-7.5 x 10~ 3 


2 


-5.859951 [2 x 10~ 5 ] 





-5.859925 [2 x 10~ 5 ] 





3 


-5.860397 [8 x 10" 6 ] 


-7.6 x 10~ 5 


-5.860450 [6 x 10~ 6 ] 


-9.0 x 10~ 5 


4 


-5.860430 [7 x 10" 5 ] 


-8.2 x 10~ 5 


-5.860380 [5 x 10" 5 ] 


-7.8 x 10~ 5 



TABLE I: Data demonstrating how sensitive the estimation of the large-/ contribution is to the number of terms N included in 
the fitting formula (|45|) . We display here results for the strong-field case (ro = 6M) and for the weak-field case (ro = 100M), 
and in both cases use the six numerically-computed modes 10 < I < 15 to fit for the part of the mode-sum with I > 15. 
[F[ >15 ]~ and [-/ ? i ' >1 5]+ are the estimated contributions from I > 15, obtained using internal and external data, respectively. 
For convenience, the values of [-F7>i5]± are given multiplied by 10 4 (M//i) 2 for r = 6M, and by 10 8 (M/» 2 for r = 100M. 
Values in square brackets indicate the fractional fitting error in [F/^sJi. The values in the 3rd and 5th columns show the 
relative differences in [-Fr>is]± with respect to the reference case N = 2, which we adopt in the rest of this work. Note that the 
differences indicated are relative to the large-/ contribution [-F/>is]± only; the differences relative to the total SF are at least 2 
orders of magnitude smaller. 

As a further robustness test for the above scheme, we check how the estimation of F[ yl would change by using 
higher multipole modes for the fitting. For this, we calculated numerically all modes up to I = 25 for ro = 6M, 100M. 
Results from this experiment are shown in Table [Til Once again we take Ff >t as the sum of all modes / > 15, but 
this time we obtain the fitting parameters Z?2n based on all sixteen modes 10 < I < 25. Again, we check how the 
results depend on N. We find that the relative difference in the value of Ff >l with respect to our reference case 
[N = 2 and fitting based on 10 < I < 15] is < 10~ 4 in all cases, as long as we take N > 2. This reassures us that it is 
sufficient to base the fitting on 10 < I < 15, as we do in the rest of this analysis. 





[*f>i 8 ]_ 


Relative difference 


[*T>vl)+ 


Relative difference 


N 


(fit using 10 < I < 25) 


w.r.t reference case 


(fit using 10 < I < 25) 


w.r.t reference case 


r = 6M 


1 


-5.106648 [1 x 10" 3 ] 


-1.2 x 10" 2 


-5.106701 [1 x 10" 3 ] 


-1.2 x 10" 2 


2 


-5.046340 [2 x 10 -5 ] 


-3.9 x 10" 5 


-5.046533 [3 x 10 -5 ] 


-6.8 x 10~ 5 


3 


-5.046634 [3 x 10~ 5 ] 


-9.7 x 10~ 5 


-5.046985 [4 x 10~ 5 ] 


-1.6 x 10" 4 


4 


-5.046500 [1 x 10" 4 ] 


-7.1 x 10" 5 


-5.047089 [1 x 10" 4 ] 


-1.8 x 10" 4 


5 


-5.046769 [4 x 10" 4 ] 


-1.2 x 10~ 4 


-5.047224 [5 x 10" 4 ] 


-2.1 x 10 -4 


r = 100M 


1 


-5.897633 [8 x 10" 4 ] 


-6.4 x 10 -3 


-5.897625 [8 x 10" 4 ] 


-6.4 x 10 -3 


2 


-5.860129 [1 x 10 -5 ] 


-3.0 x 10~ 5 


-5.860106 [1 x 10~ 5 ] 


-3.1 x 10~ 5 


3 


-5.860367 [7 x 10" 6 ] 


-7.1 x 10~ 5 


-5.860349 [8 x 10~ 6 ] 


-7.2 x 10~ 5 


4 


-5.860371 [2 x 10 -5 ] 


-7.2 x 10" 5 


-5.860304 [2 x 10 -5 ] 


-6.5 x 10~ 5 


5 


-5.860407 [9 x 10" 5 ] 


-7.8 x 10~ 5 


-5.860308 [1 x 10" 4 ] 


-6.5 x 10 -5 



TABLE II: Data demonstrating how sensitive the estimation of the high-Z contribution is to the number of modes used for the 
high-Z fitting. The structure of this table is the same as that of Table U A gain we show data for ro = 6M and ro = lOOAf , 
and [fy>i5]± describes the contribution from / > 15; however, we now use all 16 modes 10 < I < 25 for fitting. The values 
in the 3rd and 5th columns show the relative differences with respect to the reference case, i.e., N — 2 and fitting using 
10 < I < 15 — the case displayed in Table [I] This demonstrates that, taking N = 2, it is sufficient to base our large-/ fitting on 
data from 10 < I < 15. 
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Results from the above fitting procedure (with N = 2 and 10 < Z < 15) are illustrated in Fig. [4] for ro = 6M. As 
an indicative measure of the calculation error in F[ >1 we take the standard (fractional) fitting error [4l|, which we 
denote A ta ii- At least in the examples shown in Tables HI and HT1 this has roughly the same magnitude as the error 
from changing the fitting method, so A ta ii provides a realistic estimate of the total (fractional) error in determining 
the tail contribution F[ >t . The relative error from determining the large-/ tail, expressed as a fraction of the total 
SF, is 



*tail,rcl 



itail 



F, 



l>l„ 



Kin 



(48) 



In practice we find (see Appendix[DJ Ff >t /F l 



Kl„ 



io- 



T0 4 for ro in the range 6M-150M. The relative fitting 



error A ta ii,rci is a mere ~ 7 x 10 5 at tq = 6M, dropping down to ~ 2 x 10 7 at ro = 150M. 
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FIG. 4: Analytic fitting for the large-Z tail of the SF, exemplified here for ro = 6M. We used the fitting formula (I45|l with 
N — 2, and based on the modes 10 < I < 15. Circles ('©') represent actual data obtained for F^ g (calculated from r^), for the 
various modes 3 < I < 25. The dashed line is the analytic fit. Left/right panels show the same data on linear/log scales. The 
large-Z tail of the mode-sum series shows the l~ 2 fall-off expected from theory (cf. Fig. [7]). 



IV. CODE VALIDATION 



A few validation tests for our metric perturbation code were presented in BL. These included (i) a demonstration 
that the numerical solutions for the various /jM im ' s converge quadratically as h — > 0; (ii) a demonstration that these 
solutions satisfy the Lorenz gauge conditions; (iii) a confirmation that the flux of energy radiated to null infinity in 
the various modes, as calculated from our Lorenz gauge solutions, compares well with the flux obtained using other 
methods/gauges. Here we present some more validation tests, focusing on the new ingredient of the analysis, i.e., 
the calculation of the SF. We will demonstrate (i) quadratic numerical convergence of the computed SF; (ii) that the 
SF does not depend on our choice of initial data; (iii) that the full-force modes have large-Z behavior as predicted in 
theory [Eq. (|2H)) ] ; (iv) that the two one-sided values obtained for the final SF (from r$ and from r$ ) agree; and (v) 
that the total flux of energy to infinity and through the horizon is consistent with the value obtained for the temporal 
component of the SF. 



A. Numerical convergence 

The scheme introduced in Sec. IIIII should yield the final SF with a numerical error scaling as ~ h 2 (where h is the 
step size in both u and v). To check this, we performed the following test, for a selection of ro values in the range 
6M-150M and I values in the range 0-15. For given r and I we calculated the regularized mode Ffj g through the 
scheme described in Sec. IIIII (including the extrapolation to h — * 0). We recorded the values of the force calculated 
with the different resolutions hi [see Eq. Denoting these by F™ [hi], we then plotted the difference F^ g [hi]—F^ g 

(where the second term is the extrapolated force) as a function of hi- Figure [5] shows the two modes / = 2,15 at 
ro = 6 M. In each case we plot both one-sided values of the difference F^J [hi] — F^ . In all cases we find that 
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this difference decreases approximateiy like h 2 at small h, demonstrating quadratic numerical convergence. Similar 
convergence is observed for the t component. 




0.01 0.1 1 0.01 0.1 1 

h [in units of M] h [in units of M] 

FIG. 5: Numerical convergence of the calculated SF, demonstrated here for the r component, for ro = 6M. The left and right 
panels show I — 2 and I = 15, respectively. Plotted is the difference F^ g [hi] — F^„ between the value of the regularized mode 
computed with step size hi, and the value extrapolated to h — > 0. Each panel displays both one-sided values of the force: "left" 
and "right" stand for and rj values, respectively. The reference lines (dotted) have slops oc h 2 . This demonstrates the 
quadratic convergence of the numerical calculation. 



B. Dependence on initial conditions 

As explained above, we start the evolution of each of the modes h^ lm with null values along the initial surfaces, 
v = Vq and u = uq. This creates a 'spark' of spurious radiation which propagates through the grid, but dies off at late 
time. During the early transient period the numerical solution is not stationary; As the spurious waves die off, the 
numerical solution approaches its physical, stationary value. This behavior is demonstrated in the plots in Fig. (TJ5]). 
The SF has to be calculated at late enough time, to assure that the error due to residual initial waves is negligible. 
This sets a lower limit on the required evolution time T cvo , which in practice depends on ro: Waves from an initial 
spark at ro = 100M dissipate more efficiently (faster) than waves from an initial spark at ro = 6M, since the former 
experience less scattering off spacetime curvature. 

In our analysis we determined T cvo experimentally, for each value of ro considered. To assess the effect of residual 
waves, we compared the values obtained for the final SF at two different evolution times, T cvo and 0.8 x T cvo . We 
regard this difference as indicative of the error from non-stationarity in our calculation. The errors calculated for the 
various values of ro are shown in Tables HVl and fVl of Sec. fVl below. The fractional error is smaller than 10 -4 in all 
cases, hence smaller than our standard discretization error. Table Hill lists the evolution times T evo (ro) used in our 
analysis. 



C. Large / behaviour of the full- force modes 

The fact that our numerically-derived full modes [-F^'u^o^i exhibit the right behaviour at large I, through three 
leading terms in the 1/L expansion [O(L) to 0(1/ L)] provides a very strong quantitative check on our code. The 
plots in Fig. [7] demonstrate that the regularized modes F^} g = [Ff^ n (xo)]± — A±L — B a fall off faster than 1/1 at large 
/: The r component falls off as ~ L~ 2 , and the t component falls off exponentially. For the r component this indicates 
that the "singular" (large I) part of the calculated full modes is correctly described by the analytic regularization 
function A r ± L + B r , through 0(1/L). Of course, this agreement is necessary for a successful implementation of the 
mode-sum formula: If the modes [i 7 'f' u \i(xo)]± were inconsistent with the regularization function, the sum over modes 
would diverge. Note also that the procedure for evaluating the large I tail (Sec. IIIIEj) relies on the computed modes 
having the correct behaviour at large I. 
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time t (in T orbit ) time t (in T orbit ) 

FIG. 6: Time evolution of the metric perturbation. We plot the absolute values of the various functions h^ lm , evaluated at 
the location of the particle, as a function of t, for I — 2 and m = 1, 2. Top two figures are for ro = 6M; lower two figures are 
for ro = 100M. The time t is indicated along the horizontal axis in units of the orbital period, Tab- The initial transient phase 
is due to the imperfection of the initial conditions; these spurious waves dissipate rapidly, clearing the stage for the physical, 
stationary solution. The SF is calculated at late time, when the effect of the initial spurious waves is negligible. 



D. Agreement between one-sided values of the SF 

In this analysis we applied the standard version of the mode-sum formula, Eq. (|20p . rather than the "averaged" 
version, Eq. (|23p . For the r component we carried out two independent calculations of the SF, once using the internal 
values [F[^ n (xo)]-, and again using the external values [FfJ Lll (xo)]+ . (The two one-sided values of the t component 
coincide automatically, as the full modes of this component are continuous at the worldline.) This allows for an 
important consistency check: The external and internal values of the final SF must agree. Any difference between the 
two values is due to numerical error. Denoting the computational difference between external and internal values of 
the final SF by A±, we find experimentally 

A±/F r ~ i(r 5 -icr 9 (49) 

(depending on ro). This is well under the numerical error from discretisation. The experimental values of A±, for the 
different orbital radii considered, can be found in Table fVl of Sec. fVl 



E. Energy balance 

Equation relates the temporal component of the SF to the momentary rate of change of the (specific) orbital 
energy parameter £ . In terms of time t, this relation becomes £ = —(^,u t )~ 1 F t , where an overdot denotes d/dt, and 
Uq = (1 — 3A//ro) -1 / 2 . If we assume that fi/M is small enough that radiation reaction is negligible over an orbital 
period T or b ("adiabatic approximation"), then, for a circular orbit, £ also represents the average rate of change of 
£ over T or b. This must be balanced by the flux of gravitational- wave energy radiated to infinity and through the 
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r /M 


^evo/^orb 


T evo 1 ?orb 




(r component) 


(t component) 


6-10 


3 


3 


11-12 


2.5 


2.8 


13-14 


2 


2.8 


15 


1.5 


2.5 


20 


1 


2 


30 


0.8 


1.8 


40 


0.6 


1.7 


50 


0.5 


1.5 


60 


0.45 


1.5 


70 


0.4 


1.2 


80 


0.3 


1.0 


90 


0.25 


1.0 


100 


0.2 


0.8 


120 


0.15 


0.8 


150 


0.12 


0.6 



TABLE III: Evolution times taken in our analysis. These were chosen (experimentally) long enough to guarantee that any 
residual effect from the initial spurious waves is negligible. T cvo is the numerical evolution time, and T or b = 2ir/Qo is the orbital 
period. 
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FIG. 7: Large-/ behavior of the regularized modes F£L g . The reference lines in the left panel are oc I 2 ; the reference line in 
the right panel is oc exp( — 1.22/). The left panel displays the radial component F^} s (calculated using internal derivatives) for 
3 < I < 25 and for ro = 6M, 100M. The regularized modes appear to fall off as oc l~ 2 at large I, in agreement with theory. 
The right panel shows the /-modes of the t component (for ro = 6M); these are already "regularized", and show an exponential 
fall-off at large /. (Note the scale for the two panels is different: log-log for the left panel, semi-log for the right.) 



horizon, averaged over T or b. If we denote the former by Eqo and the latter by -Eeh (both taken positive), we have the 
energy balance formula 

Stotal = Eoo + E E r = -[it = Ft/ul (50) 

Both asymptotic fluxes and E^r can be constructed from the perturbation fields h^ Lm , evaluated at the corre- 
sponding asymptotic domains. Validity of Eq. (|50|) then provides a strong qualitative test of our calculation. 

We can readily express the asymptotic fluxes E x and -Eeh in terms of the h^' lm, s, with the help of the Weyl 
scalars, = ~C a f3jsl a m^Pm 5 and V>4 = —C a ^sn a m* ,3 n l m* s . Here C Q/ 3 7 «5 is the Weyl tensor corresponding to the 
perturbation h a p, and l a , n a , and m a are the Kinnersley null vectors, given by l a = 1, 0, 0) n a = i(l, — /, 0, 0), 
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and m a — -^:(0, 0, 1, jj^p). Decomposing the perturbation as in Eq. ([16]) . we obtain the asymptotic relations 



oo I 

M r °°) = E E 

1=2 m=-l 
oo I 

Mr - 2M) = E E 

Z=2 m=-i 



1 



4Z(Z + l)Ar 
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i ft 



(8) 



,(9) 



f(7) -(10) 
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{D 2 Y lm + i{sm6)- x D 1 Y l 



lm ) 



(51) 



The first relation is valid at leading order in 1/r, and the second at leading order in /. In obtaining these relations 
we have made the replacements fd r — > —dt (for ip^ at null infinity) and fd r — > dt (for tpo at the horizon). For circular 
orbits, the asymptotic fluxes are given in terms of the Weyl scalars as [HI, |43| 
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(53) 



where the integration is carried out over 2-spheres r = const — » oo and r = 2M, respectively. Now proceed as follows: 
(i) Substitute Eqs. ([ST]) in Eqs. (|55]l. (ii) For £eh apply the asymptotic gauge conditions ft' 2 ) = ftW, ft' 4 ) = ftA% and 
ft( 8 ) = ft( 9 ) [see Eqs. (|A13[I - (|A16|) in Appendix [A"] . (iii) Replace d/dt — > —imflo. (iv) Integrate over the spheres using 
the formulas (A4) of BL. This yields the final relations 
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and ftgjj are these fields 



where A = (i — 1)(Z + 2), ft^ are the fields ftW im evaluated at null infinity (u <C v — > oo) orir:| 
evaluated at the event horizon (v <C u — > oo). Eq. ([54]) for _Eoo agrees with Eq. (57) of BL, which was derived directly 
from the Issacson effective energy-momentum tensor. 

To test our calculation of F t , we used our evolution code to obtain the energy fluxes at infinity and through the 
horizon, based on Eqs. ([54]) and ([55]) . and then checked consistency with the balance equation (150|) . To extract we 
evaluated the numerical solutions ft( 7 < 10 ) at v = 5200M, 9000M, 15000M and u = 800M, 3000M, 5000M, for orbital 
radii in the ranges 6M < tq < 20M, 20M < ro < 100M, and tq > lOOAf, respectively. To derive £eh we evaluated 
^(1,5,7,9,10) a t u = 5200M, 6500M, 7500M and v = 800Af, 1500M, 2500M for the above corresponding values of r . 
These values were selected experimentally such that the fractional error in the total flux (from the finite extraction 
distance and from the spurious initial waves) is less than 10 -4 for each of the / modes. We then used these values in 
Eqs. ([54"]) and ([55]) , summing from I = 2 to I = i ma x, where Z max was determined experimentally, requiring that the 
fractional truncation error in the total flux (from omitting the modes / > l max ) is < 10~ 4 . This required Z max values 
between 9 (for r = 6 M) and 4 (for r = 150M). 

Following the above procedure, we obtained and -Eeh for a list of orbital radii between 6M and 150M. The 
ratio -Eeh/^oo turns out very small for all radii, decreasing monotonically with ro from 3.3 x 10~ 3 at ro = 6Af to 
150M (these values are consistent with Martel's [40]). The values obtained for the total energy flux, 



-Etotai = Eoo + Eeh, are listed in Table ITVl of Sec. IVl below. We find that the fractional difference 
less than ~ 5 x 10 -4 in all cases, providing a strong quantitative check of our results. 



jEtot&i / F t — 1 



is 



20 



r /M 




(M/fi) 2 F t /u 


t 





rel. 


diff. 


6.0 


-1.99476 


X 


10" 


-3 


[7 


X 


1(T 


5| 


[6 


X 


10" 


61 


9.40338 


X 


10 


-4 


9.40190 


x 


10" 


-4 


1 


6 


X 


10" 


4 


6.2 


-1.60515 


X 


10" 


-3 


[8 


X 


10" 


5] 


[1 

L 


X 


10" 


Si 


7.81183 


X 


10 


-4 


7.81064 


x 


10" 


-4 


1 


5 


X 


10" 


4 


6.4 


-1.30550 


X 


10" 


-3 


[7 


X 


10" 


5l 


[1 

L 


X 


10" 


5l 


6.54180 


X 


10 


-4 


6.54101 


X 


10" 


-4 


1 


2 


X 


10" 


4 


6.6 


-1.07197 


X 


10" 


-3 


[7 


X 


10" 


5i 


[1 

L 


X 


10" 


Si 


5.51794 


X 


10 


-4 


5.51723 


X 


10" 


-4 


1 


3 


X 


10" 


4 


6.8 


-8.87844 


X 


10" 


-4 


[6 


X 


10" 


5] 


[1 


X 


10" 


Si 


4.68497 


X 


10 


-4 


4.68411 


X 


10" 


-4 


1 


8 


X 


10" 


4 


7.0 


-7.41101 


X 


10" 


-4 


[7 


X 


10" 


5l 


[1 


X 


10" 


Si 


4.00157 


X 


10 


-4 


4.00117 


X 


10" 


-4 


1 





X 


10" 


4 


7.2 


-6.23065 


X 


10" 


-4 


[7 


X 


10" 


5i 


[1 

L 


X 


10" 


Si 


3.43687 


X 


10 


-4 


3.43627 


X 


10" 


-4 


1 


7 


X 


10" 


4 


7.4 


-5.27271 


X 


10" 


-4 


[6 


X 


10" 


5 1 


[1 


X 


10" 


Si 


2.96692 


X 


10 


-4 


2.96645 


X 


10" 


-4 


1 


6 


X 


10" 


4 


7.6 


-4.48905 


X 


10" 


-4 


[6 


X 


10" 


5] 


[1 


X 


10" 


Sl 


2.57336 


X 


10 


-4 


2.57288 


X 


10" 


-4 


1 


9 


X 


10" 


4 


7.8 


-3.84324 


X 


10" 


-4 


[5 


X 


10" 


5i 


[1 

L 


X 


10" 


Si 


2.24184 


X 


10 


-4 


2.24148 


X 


10" 


-4 


1 


6 


X 


10" 


4 


8.0 


-3.30740 


X 


10" 


-4 


[5 


X 


10" 


5] 


[1 

L 


X 


10" 


Sl 


1.96105 


X 


10 


-4 


1.96066 


X 


10" 


-4 


2 





X 


10" 


4 


9.0 


-1.66810 


X 


10" 


-4 


[5 


X 


10" 


5i 


[1 


X 


10" 


Sl 


1.05933 


X 


10 


-4 


1.05908 


X 


10" 


-4 


2 


4 


X 


10" 


4 


10.0 


-9.19067 


X 


10" 


-5 


[3 


X 


10" 


5] 


[9 


X 


10" 


61 


6.15158 


X 


10 


-5 


6.15047 


X 


10" 


-5 


1 


8 


X 


10" 


4 


11.0 


-5.41605 


X 


10" 


-5 


[3 


X 


10" 


5] 


[2 

L 


X 


10" 


Sl 


3.77904 


X 


10 


-5 


3.77856 


X 


10" 


-5 


1 


3 


X 


10" 


4 


12.0 


-3.36587 


X 


10" 


-5 


[2 


X 


10" 


5] 


[2 

L 


X 


10" 


Sl 


2.42911 


X 


10 


-5 


2.42857 


X 


10" 


-5 


2 


2 


X 


10" 


4 


13.0 


-2.18388 


X 


10" 


-5 


[2 


X 


10" 


5] 


[2 


X 


10" 


Sl 


1.62071 


X 


10 


-5 


1.62022 


X 


10" 


-5 


3 


1 


X 


10" 


4 


14.0 


-1.46851 


X 


10" 


-5 


[1 


X 


10" 


4i 


L 


X 


10" 


Sl 


1.11574 


X 


10 


-5 


1.11564 


X 


10" 


-5 


8 


5 


X 


10" 


5 


15.0 


-1.01772 


X 


10" 


-5 


[1 


X 


10" 


4i 


f2 

L 


X 


10" 


Sl 


7.88904 


X 


10 


-6 


7.88597 


X 


10" 


-6 


3 


9 


X 


10" 


4 


20.0 


-2.25549 


X 


10" 


-6 


[6 


X 


10" 


5"i 


[6 


X 


10" 


61 


1.87151 


X 


10 


-6 


1.87111 


X 


10" 


-6 


2 


2 


X 


10" 


4 


30.0 


-2.80813 


X 


10" 


-7 


[4 


X 


10" 


5 1 


[4 


X 


10" 


Sl 


2.48643 


X 


10 


-7 


2.48600 


X 


10" 


-7 


1 


7 


X 


10" 


4 


40.0 


-6.51219 


X 


10" 


-8 


[3 


X 


10" 


5i 


L 


X 


10" 


61 


5.95007 


X 


10 


-8 


5.94897 


X 


10" 


-8 


1 


8 


X 


10" 


4 


50.0 


-2.10849 


X 


10" 


-8 


[2 


X 


10" 


5"i 


[4 


X 


10" 


Sl 


1.96249 


X 


10 


-8 


1.96203 


X 


10" 


-8 


2 


3 


X 


10" 


4 


60.0 


-8.41306 


X 


10" 


-9 


[9 


X 


10" 


5] 


[3 


X 


10" 


Sl 


7.92670 


X 


10 


-9 


7.92424 


X 


10" 


-9 


3 


1 


X 


10" 


4 


70.0 


-3.87411 


X 


10" 


-9 


[8 


X 


10" 


5] 


[4 


X 


10" 


5l 


3.68189 


X 


10 


-9 


3.68086 


X 


10" 


-9 


2 


8 


X 


10" 


4 


80.0 


-1.98069 


X 


10" 


-9 


[7 


X 


10" 


5"i 


[S 


X 


10" 


Sl 


1.89462 


X 


10 


-9 


1.89360 


X 


10" 


-9 


5 


4 


X 


10" 


4 


90.0 


-1.09654 


X 


10" 


-9 


[6 


X 


10" 


5] 


[6 


X 


10" 


Sl 


1.05415 


X 


10 


-9 


1.05365 


X 


10" 


-9 


4 


8 


X 


10" 


4 


100.0 


-6.46305 x 10" 


10 


[6 x 10" 


- 5 ] 


[4 x 10" 


- s ] 


6.23806 x 10" 


10 


6.23628 x 10" 


10 


2 


9 


X 


10" 


4 


120.0 


-2.59096 x 10" 


10 


[5 x 10" 


- 5 ] 


[3 x 10" 


- s ] 


2.51573 x 10" 


10 


2.51496 x 10" 


10 


3 


1 


X 


10" 


4 


150.0 


-8.47172 x 10" 


11 


[4 x 10" 


- 5 ] 


[6 x 10" 


- 8 ] 


8.27475 x 10" 


11 


8.27279 x 10" 


11 


2 


4 


X 


10" 


4 



TABLE IV: The temporal component of the SF, as a function of the orbital radius ro. Values in the first square brackets in 
the second column are estimates of the fractional numerical error in F l from the finite-grid discretization, A*[ iscr [see Eq. (|42[) ], 
Values in the second square brackets in the second column are estimates of the fractional error from residual non-stationarity 
of the late-time numerical evolution (which is mainly due to leftover spurious waves arising from the imperfect initial data). 
The third and fourth columns compare between the work done by the temporal SF and the total flux of energy radiated in 
gravitational waves, the latter extracted from the numerical solutions using the procedure described in Sec. IIV El The last 

column displays the relative difference |i5totai/(Fi/iio) — lj, showing that the balance equation (|50[) is satisfied within the 

numerical accuracy, and providing a strong quantitative check of our results. 



V. RESULTS 



A. Temporal component 

We calculated the t component of the SF for 29 values of the orbital radius, in the range from ro = 6M to 
ro = 150M, using the procedure described in Sec. I1III The results are displayed in Table HVl The computation error 
in F l is estimated at < 10~ 4 for all radii considered. The table also shows, for each of the radii considered, how the 
work done by the temporal component of the local SF is balanced by the total flux of radiated energy. 
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B. Radial component 

We calculated the radial component of the SF for 29 values of the orbital radius, in the range from ro = 6M 
to to = 150M. Tables IVII and IVIII in Appendix [D] present results for the internal and external values of the SF, 
respectively (recall the two one-sided values are expected to agree with each other, within numerical error). In each 
table we indicate separately the two contributions [F[ <15 ]± and [F[ >15 }±, where the former is the part calculated 
directly using our evolution code, and the latter is the extrapolated contribution from I > 15, calculated as explained 
in Sec. IIII El above. The large-l tail contributes at most 2% of the total SF (depending on ro). The relative fractional 
error in [-F;> 15 ]± [he., A ta ii,roi, calculated through Eq. (|48p ] is at most comparable to (and mostly much smaller than) 
the fractional discretization error in [F[ <15 }±, which itself is at most ~ 10~ 3 . 

In Table IVl we present our final results for the radial component: the SF as a function of the orbital radius ro. As 
the 'final' result we quote the average between the two (nearly identical) one-sided values. The third column displays 
the magnitude of the difference between the two one-sided values, which is entirely due to computational error. The 
relative magnitude of this error (given in square brackets in the third column) is in all cases much smaller than the 
fractional discretization error in F r . The latter error, given in square bracket in the second column, is taken as the 
average between the two one-sided discterization errors Aj iscr estimated from Eq. (j42| . 

The last column of Table [V] displays the estimated error from residual non-stationarity of the numerical solutions 
(from residual spurious initial waves). Displayed is the difference in the values F r obtained at different evolution times 
(as explained in more detail in Sec. IIV 51 above) . The values in square brackets in the last column show the fractional 
error from non-stationarity relative to the total SF. This error is in all cases much smaller than the discretization 
error. 

Thus, the dominant source of error in our analysis is associated with the finite-grid discretization of the field 
equations. We estimate our final results for F r to be correct to within at least ~ 0.1% for all orbital radii considered. 
The results for 10M < ro < 20M are likely to be correct to within ~ 0.01%, and the results for ro > 20M to within 
mere ~ 0.001%. 

We plot F r (ro) in Fig. [8l The radial SF is "repulsive" (i.e, acting outward, away from the central black hole) for 
all ro. At large orbital radii the numerical data can be fitted analytically as 



F r (r » M) ~ ^ 



M 

ao + oi — 

7"0 



M 
a 2 1 — 

ro 



M 
a>3 | — 



with 



a = 1.999991, a 1 = -6.9969, a 2 = 6.29, a 3 



-24.6. 



(56) 



(57) 



This formula reproduces the numerical data within the numerical accuracy [< 10 -3 ] for all ro > 8M. The leading- 
order term, F r ~ ao/x 2 jr\ ~ 2/x 2 /ro is consistent with the "Keplerian" SF describing the back-reaction effect from the 
motion of the black hole about the system's center of mass (we discuss this below, when considering the SF correction 
to the orbital frequency). Near the Innermost Stable Circular Orbit (ISCO), ro = 6M, we fit the numerical data 
analytically as 



F r (r > 6M) ~ ^ (1 - 2M/r ) (b + hx + b 2 x 2 Q + b 3 x 3 ) , 



(58) 



where xo = 1 — QM /ro and the 'best fit' parameters (based on data in 6M < ro < 8M) are given by 

b = 1.32120, bi = 1.2391, b 2 = -1.297, b 3 = 1.07. 
This reproduces the numerical data within the numerical accuracy for all ro < 8M.(49| 



(59) 



C. Conservative shift in the orbital parameters 



Given F r , we can calculate the shift in the orbital energy and angular momentum parameters using Eq. ([5]). The 
relative shifts AS = (S—So)/Sq and AC = (£— Co)/ Co are plotted in Fig.0 Recall that this effect is gauge-dependent; 
the values computed here for AS and AC are the Lorenz-gauge values. 

The shift in the orbital frequency can be derived from Eq. ([5]). At large ro we obtain, using Eq. (|56p . 
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TABLE V: Final results for the radial component of the SF. The second column lists the values obtained for the various orbital 
radii ro, taken as the average between internal and external values: ([-F r ]+ + [F r ]_)/2. Values in square brackets in the second 
column are estimates of the fractional numerical error in F r from the finite-grid discretization (see the text for details) . The 
third and fourth columns display estimates of the magnitude of error from two other sources: The third column shows the 
magnitude of the difference — [F 7 ]_, which is entirely due to numerical error; the values in square brackets give the 

fractional error 2([F r ]+ — [F r ]-)/([F r ]+ + [F r ]-). The fourth column shows the estimated error from residual non-stationarity 
of the late-time numerical evolution, which is mainly due to leftover spurious waves arising from the imperfect initial data; 
values in square brackets again describe the fractional error. Both sources of error contribute negligibly to the overall error in 
the SF, which is therefore dominated by the discretization error. 



where ci =3ao — ai,C2 = 3ai — a2, and C3 = 3a2 — 03, with the coefficients a n given in Eq. (|57[) . The term proportional 
to ao (~ 2) is precisely the "Newtonian" SF [see, e.g., Eq. (2) of 35]], which dominates the SF effect at ro 3> M. 
This piece of the force is simply the 0(/u) difference between the standard Keplerian frequency fl 2 = (M + n)/R 3 
(expressed in terms of the separation R) and SIq = M/tq, with the separation R related to the "center-of-mass" 
distance ro through M(R — ro) = fj.ro. The rest of the terms in Eq. (|60[) are general relativistic (3PN) corrections. 
We define A(0 2 ) G r = (ft 2 - ftg)/fto + W M > and in Fig.[UJ]plot A(ft 2 ) GR as a function of r . 

It is natural to ask how our result for ft compares, at large ro, with results from PN literature. We must first 
note that, although ft is gauge- invariant, r itself is not, and so the functional form of ft(ro) in Eq. ([60]) is gauge 
dependent. This makes it difficult to compare with the standard (non-perturbative) PN result [44]], which is given 
in a particular coordinate gauge (the "harmonic" gauge) that, in a perturbative context, does not coincide with the 
Lorenz gauge employed here. Another calculation of the conservative PN SF was carried out recently by Nakano [4|| , 
within perturbation theory, using a modified version of the Regge-Wheelar gauge. The results from this calculation, 
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FIG. 8: The radial component of the SF. The data in the upper panel correspond to the second column in Table IVl The inset 
shows an expansion of the ISCO area. The solid line is a plot of the large-ro analytic fit given in Eq. (|56p . The dashed line 
represents the complementary near-ISCO fit (158ft . The lower panel shows the relative difference between the numerical data and 
values obtained using the analytic fit formulas. Using the large-ro fit for ro > 8M and the near-ISCO fit for 6M < ro < 8M, 
one recovers all numerical data to within the numerical accuracy. 




FIG. 9: The shift in the energy (left panel) and angular momentum (right panel) parameters, due to the conservative SF effect. 
The plots show the relative shifts AS = {£ — £q)/£q and AC = (£ — Co) /Co, where £o and Co are the geodesic values. Solid 
and dashed lines correspond to the large-ro and near-ISCO analytic fits, Eqs. (I56|l and (|58|) . respectively. Recall AS and AC 
are gauge dependent; the values shown here are the Lorenz-gauge values. 



too, cannot be directly compared to ours, because of the different gauges used. 

One may attempt to circumvent the gauge ambiguity problem by writing down an expression (in a PN form) for one 
gauge-invariant quantity in terms of a second gauge-invariant quantity — such an expression would be 'truly' gauge 
invariant and would allow direct comparison between calculations done in different gauges. For circular orbits, both fi 
and S = £ — flC are gauge invariant (see Sec. Ill B]) . and we may attempt an expression of the form S(Q). Introducing 
the new gauge-invariant variable x = (Mil) 1 ' 3 , we obtain, at 3PN, 

S = 1 - §^ - fa* - fg* a + 0(/* a ), (61) 

with a vanishing 0(/i) term. Hence, S(x) is not useful for comparing the conservative SF effect in the case of a 
circular orbit. Comparison of our results with results from PN literature is not at all straightforward, and we leave it 
for future work. 
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FIG. 10: Shift in the orbital frequency due to the conservative SF effect. Shown is the quantity A(fi 2 )GR = (fi 2 — nj^/n 2 , — 
( — 2/i/M), describing the relative shift in Q. 2 , minus the leading-order Keplerian term. The solid line shows the large-ro (3PN) 
analytic fit, Eq. (|60[) . The dashed line corresponds to the near-ISCO fit (|58[) . 



VI. CONCLUDING REMARKS AND FUTURE APPLICATIONS 



This work marks a minor milestone in a long-term program aimed to develop the theoretical and practical tools 
for computing EMRI orbits (and, eventually, their gravitational waveforms). We compute here for the first time the 
gravitational SF in an example of a particle orbiting a black hole, demonstrating the applicability of our approach, 
whose main elements are (i) direct solution for the metric perturbation, in the Lorenz gauge; (ii) numerical evolution 
in the time domain; and (hi) use of the mode-sum scheme to derive the local SF. In the case of a strictly circular orbit, 
the analysis of the local SF provides us with little new physics: The radiative effect is well known from energy-balance 
analysis, and the conservative force does not have a strict gauge-invariant significance. Calculation of gauge invariant 
conservative effects (like the shift in the ISCO frequency, or the correction to the rate of perihelion precession) requires 
analysis of (at least slightly) non-circular orbits. In follow-up work we intend to extend our analysis to eccentric orbits 
(see below), which would gain us access to this more interesting physics. 

Self-force calculations bring about major issues of computational cost and computational efficiency. All compu- 
tations in this work were carried out on a standard desktop computer (3GHz dual- processor, with 4Gb of RAM). 
Calculation of the SF at a single strong-held point, with fractional accuracy < 10~ 3 , took ~ 2 hours of CPU time. 
This is practical enough for studying the simple one-parameter family of circular orbits, but may not be practical for 
studying more general orbits. There are a few obvious ways by which one may improve the efficiency of the numerical 
algorithm: (1) Our evolution code currently utilizes a uniform grid. This is very inefficient, since the resolution 
requirement near the worldlinc is much higher than anywhere else on the 2-dimensional grid. Our problem naturally 
calls for a mesh- refinement treatment. This is a standard technique in numerical relativity, but its application would 
require a major modification of our code. (2) We may try to improve the rate of decay of the initial spurious waves, by 
using the stationary numerical solutions obtained with low resolution as initial conditions for the evolution at higher 
resolution. This will allow to evolve for shorter periods, hence saving computation time. It should be straightforward 
to implement such a procedure. (3) In the present analysis we conservatively set the same accuracy threshold for 
each individual I mode of the force. Since the contribution of the individual modes to the total SF vary over a few 
orders of magnitude, this procedure is not very economic. It would be better to use an algorithm which incorporates 
a threshold on the total force. This, too, could be implemented rather easily. 

Since our code is based on time-domain evolution (with no frequency decomposition), it is readily extensible to deal 
with any orbit in Schwarzschild spacetime. The finite-difference algorithm would change slightly (it would resemble 
the algorithm used for radial plunge trajectories [H, H(|), but the stability features and resolution requirements of 
the code would not change. Work to extend our analysis to eccentric orbits is now in progress. 

It is more challenging to apply our approach for orbits in Kerr spacetime. In this case we may no longer rely on 
a spherical-harmonic decomposition of the held equations, and — insisting on a time-domain analysis in the Lorenz 
gauge — we would have to apply time evolution in 2+1-D. The challenge here is two-fold: Firstly, the solutions to the 
2+1-D held equations are no longer continuous along the worldline (as in the 1+1-D case), but rather diverge there 
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logarithmically. Secondly, a stable numerical scheme for evolution of Lorenz-gauge perturbations in 2 + ID is yet to 
be developed. A numerical scheme for dealing with the first of the above difficulties had been outlined in Sec. V of 
BL, and was recently implemented for a scalar- field toy model (46l |. 
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APPENDIX A: FIELD EQUATIONS AND GAUGE CONDITIONS FOR THE PERTURBATION 

FUNCTIONS h (i)lm (r,t) 

We give here explicit expressions for the various terms appearing in our basic set of mode-decomposed field equations 
([19]). We use the notation / = 1 - 2M/r, f = 2M/r 2 , f = l- 2M/r , and X — (I + 2) (I - 1). d r is taken with 
fixed t, and d v is taken with fixed u (v and u are the standard Eddington-Finkelstein coordinates). For brevity we 
occasionally omit the indices I, m. 

The terms M ^h^ in Eqs. (| L9[) are given by 



M 



M%m = \ffhf + ^(1 - 4M/r) (fcM - ^ 5 > - fh&) - - 6M/r)h^, 

- h u) = l _ ff - h ( ? + r _ h nj + f_ ^( 2) _ m ^ _ fj_ _ m _ 03) _ 2fU ^ 
M { %m = —L \m - m - (i - 4M/o ( + m 



M%h& = \r (h$ - nf) + 1) (//r 2 )^ 2 ) - \ff/r [zm + 2m - m + in + ^ 

M %^ U) = 4 ( X - 4.5M/r)^ (5) - + 1) (h {1] - A (3) ) + hi - 3M/r) (z(/ + l)h {6] - /i (7) ) 



_ /jC5) _ (! _ 4M / r ) ^0) + h (6)\ 

= If (h§) - b$) - \ff/r (Zh^ + 2m - h m ) , 
M ( f 3) m = 4 (1 - 4.5M/r) - JL (1 - 3M/r) ^ 10 \ 



(8) 
(j) 



2r 2 



A1 ( ( 10)^) = _i_^(10) +A ^(9)j. 



(Al) 
(A2) 
(A3) 
(A4) 
(A5) 
(A6) 
(A7) 
(A8) 
(A9) 
(A10) 



For a circular equatorial geodesic orbit with r = r$ [hence with i-frequency Qo — (M/r^) 1 / 2 and specific energy 
£q = (1- 2M/r )(l - iM/ra)- 1 / 2 ], the source terms S^ lm in Eqs. JTSJ) read 



S(0i"»( rj t) = 47r5 a (i) <5(r-r ) x 



r im *(7r/2,r2 i): * = 1—7 (even parity modes), 
Y l e m *(ir/2, Sl t), i = 8—10 (odd parity modes). 



(All) 
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The coefficients are given by 

- / 2 Ao, 
« (3) = /o/ro, 

a (2) = a (5) = a (9) = 0i 

« (6) = r nl 

a ( 7 ) = r n 2 [l(l + 1) - 2m 2 ], 

a {10) = 2imr Q n 2 l . (A12) 

The four Lorenz gauge conditions h a 0^ — translate, upon decomposing in tensor harmonics, to four constraints 
on the time-radial functions ftw. These read 

- hf + f (-hf + hf + m/r - m/r) = 0, (A13) 

hf - 0$ + fhf - (f/r) (h^ - _ /^(3) _ 2 /fcW) = 0, (A14) 

/ (hf + 2h {5) /r + 1(1 + 1) h {6) /r - h {7) /A = 0, (A15) 

/ (h^ + 2h {9) It - h {1 V /r) = 0. (A16) 



(4) 



APPENDIX B: FORMULAS FOR THE COEFFICIENTS f n± AND f n± 



We give here formulas for constructing the various coefficients /" ± appearing in Eq. (J29J), using the Lorenz-gauge 
metric perturbation fields h^ lm (r,t) and their derivatives. For brevity we omit the superscripts l,m from both the 
fn± l, s and the h^ lm, s. We use here the notation C = A)Aoj where, recall, Cq is given in Eq. ((4]). Also recall So 
is given in Eq. ((3]), /o = (1 — 2M/ro)_, and X = (I + 2)(l — 1). r derivatives are taken with fixed t, and t derivatives 
are taken with fixed r. All functions h^ lm and their derivatives in the expressions below are evaluated at r = ro and 



t = to. Subscripts '±' refer to taking r derivatives from r 
For the r component we have 



(we omit these subscripts whenever f°, = /"_). 



1 
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(B3) 
(B4) 
(B5) 
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ft = fr = 0. (B7) 

For the t component we have 

/< = —ElLlf^hf + ±imS £ (2f - £ 2 )f^h^ + l -£ 2 C 2 f^(M/r a )h^ 

--/„ roA)(^+/o)fc,t -j»m£ fo/ ft 27(7TTj '* + 21(1 + 1) k 

imC Q £lf^(M/r ) - (5) 1 - 2 , ^ (6) 1 s.-it (6) m 2 £ 2 (£ 2 + / )/ - 2 r - (7) 

+ 4l(7TT)A 71 ' (B8) 



f t _ ^o/o 1 tf4) _ ^o( g o +/o)/o 2r O i(7) _ 5im£o^o// 7(7) 
/2 2Z(Z + 1) 41(1 + 41(1 + 1)X ' 



(BIO) 



.< £o( g o + /o>o/ 2 r(7) , iTOgo£|/g 1 - (7) 

/3 4Z(/ + 1)A kt + 41(1 + h ' (BU) 

rt im ^ofa\ (8) imCl(8 2 + /o)/ ~ 2 ro - ( i 0) m 2 ^/- 1 - (10) 

^-"WW 2/(7TT)a ^ + /(/ + i)A h (B12) 

/s= 2/(Z + l) ^ 21(1 + 1) h + 21(1 + 1)X h ' (B13) 



APPENDIX C: FORMULAS FOR THE COUPLING COEFFICIENTS 

We give here formulas for re-expanding all angular functions in Eq. ([2^)) in spherical harmonics Y lm (8,(p). The 
following identities hold for all /, m and any 9, (p. We have 

sin 2 6Y lm = a l ( ™ 2) Y l+2 > m + a l ^Y lm + a[ m 2) Y l - 2 ' m , 

cos sin 01^™ = ) Y l+2 > m + p\™ } Y lm + /?(™ ) Y l ~ 2 ' m , 

sin 2 9Y l ™ = j[^ 2) Y l+2 ^ n + j[^Y lm + j l ^ 2) Y l - 2 ' m , 

smOY 1 ^ 1 = 5|™ ) y i+1 ' m + <y[? 1) y , - 1,m , 

cos6»F' m - sm6Y l e m = e i ( ™ 1) r i+1 ' m + ^ 1) F ,_1 ' m , 

• 3 r\-\rlm Am \^l-\-3,m , /-lm \rl-\-l,7n 1 Am yi-l,m 1 /*£m — 3,m 

sin wr e - c,( +3 )X +(,(+!) j + ( >(-i) r +c »(-3) r > 

cos e sin 2 0y{$ = d™ 3) Y i+3 ' m + g™Y l+1 > m + 1) y ,_1,Tn + d™ 3 ) y'- 3 ' m . (CI) 
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The coefficients are all constructed from 



using 



lm / < f 1 lm i r~i~£ /~i~z 

a (+2) ~ ~ M+l,m^!+2,mi a (Q) — 1 — — °/+l 



r — m 2 



(2Z + 1)(2Z — f) 



In i 



1/2 



*(-2) 



-Ci m Ci—i, m , 



P(+2) — 'C;+l,mC;+2,m, /^(o) — ^Cf+l, m ~~ + l)Qmi /^(-2) — — G + l)^ImCl-l,mi 



(C2) 

(C3) 
(C4) 



7 ™ a) = re, + i, ra C 1+2 , m , 7(S) = ™ - Ki + 1) + l'Ci+i,m + + iyCf m , 7(- 2) = (I + irC| m C H , rai (C5) 



^(+1) = ^l+l,m) fy-l) — — (I + 



= (+1 



(C6) 
(C7) 



Clm 
:+3) 

/-lm 

M+i) 

/-lm 

C(_i) 

/-lm 

k-3) 



. in • 

Ci + l.m [/ (l - Cj+l, m - Cf +2 , m ) + (Z + l)C; 2 r 

— Cj,m [(' + 1) (l - Cj-1,to ~ Q 2 m) + ^J+l,rr. 
(Z + l)Ci im C;_i im C;_2,m, 



(C8) 



flm 

M+3) 

elm 

M+i) 

elm 

M-i) 

elm 
M-3) 



' 2 C| + l,mCi+2,mCl+3,m, 

c m , m [™ 2 - i(i + 1) + z 2 c 2 +l!tn + (/ + i) 2 c 2 m + z 2 c? +2 , m ; 

C i>m [™ 2 - Z(Z + f ) + Z 2 Q 2 +lim + (Z + l) 2 C 2 m + (Z + 1) 2 C 2 _ 

(Z + l) 2 C;^ m C;_i^ m C;_2,m- 



(C9) 



APPENDIX D: DETAILS OF NUMERICAL RESULTS FOR F r 

In this appendix we tabulate data obtained for the radial component F r , braking it up into low-Z and high-Z 
(extrapolated) contributions, and displaying separately internal and external values. These data is used in Sec. fVlfor 
error analysis. 
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TABLE VI: Results for the radial component of the SF, broken up into the contribution from I < 15 (computed directly using 
our numerical evolution code), and the contribution from the I > 15 tail (estimated as described in Sec. IIII E|l . Presented here 
are "internal" values of the SF (i.e., those obtained through taking one-sided derivatives of the perturbation from r — > r$). 
The "external" values (which should be the same within numerical error) are given in Table rVTll below. The values in square 
brackets in the second column represent the fractional discretization error Aji acr , as estimated using Eq. (|42|l . A ta ii,rci is the 
relative fractional error in F[ >15 (i.e., error expressed at a fraction of the total SF), which is estimated using Eq. (|48[) . In our 
analysis we made sure that the error from estimating the contribution from the large-Z tail is kept smaller than the discretization 
error. 
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TABLE VII: Same as Table [VTU here for the external values of the SF. 



Regge and J. A. Wheeler, Phys. Rev. 108, 1063 (1957). 

Mino, Prog. Theor. Phys. 99, 79 (1998). 

Mino and H. Nakano, Prog. Theor. Phys. 100, 507 (1998). 

Nakano and M. Sasaki, Prog. Theor. Phys. 105, 197 (2001) [arXiv:gr-qc/0010036 
Barack and A. Ori, Phys. Rev. D 64, 124003 (2001) [arXiv:gr-qc/0107056|. ~ 
Barack and C. O. Lousto, Phys. Rev. D 72, 104026 (2005) arXiv:gr-qc/0510019 . 

2010 . 



Detweiler and E. Poisson, Phys. Rev. D 69, 084019 (2004) |arXiv:gr-qc/0312010| . 
W. Misner, K. S. Thorne, and J. A. Wheeler, Gravitation (Freeman, San Francisco, 1973). 
Detweiler, Class. Quant. Grav. 22, S681 (2005) [arXiv:gr-qc/050l"004] . 
A. Ori, Phys. Rev. D 70, 124027 (2004) [arXiv:gr -qc/0312013|. 
F. J. Zerilli, Phys. Rev. D 2, 2141 (1970). 

C. O. Lousto and R. H. Price, Phys. Rev. D 56, 6439 (1997 ) [arXiv:gr-qc/9705071 . 
K. Martel and E. Poisson, Phys. Rev. D 66, 084001 (2002 ) |arXiv:gr-qc/0107104| . 
K. Martel, Phys. Rev. D 69, 044025 (2004) [arXiv:gr-qc/0311017| . 

W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery, Numerical Recipes (Cambridge University Press, 
Cambridge, 1992). 

S. A. Teukolsky, Astrophys. J. 185, 635 (1973). 

A. Teukolsky and W. H. Press, Astrop hys. J. 193, 443 (1974 ). 
Blanchet, Living Rev. Rel. 9, 4 (2006) |arXiv:gr-qc/02020i6]. 



Nakano, unpublished work; Presentation at the Ninth Capra meeting at UWM (2006). 
Barack and D. A. Golbourn, in preperation. 



By "multipole modes" we mean here, and throughout this paper, the components of a spacetime function in a spherical- 
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harmonic basis defined on 2-spheres centered at the black hole (and not at the particle), as usual in black hole perturbation 
theory. 

[48] We use here the same notation and definitions as in BL, except for a redefinition /i (3) [here] = (1 - 2M/r)"'/i (3> [BL]. This 
simplifies slightly the form of the perturbation equations (|19|) . and is motivated by the fact that [BLpaper] vanishes at 
the horizon. Also, we fix here an error in Eq. (20) of BL: The factor "i" in the expressions for fej,™ and h l J£ is erroneous 
and we omit it here. 

[49] The analytic-fit formula (|58[) should not be regarded as representing the first few terms in a convergent expansion in 
x = 1 — 6M/ro; In fact, it is evident from Fig. [8] (lower panel) that such a series is likely to have a very small radius of 
convergence around ro = 6M. We give the analytic fit ((58} merely for practical reasons, as it correctly gives the value of 
the SF anywhere in the strong-field regime 6M < ro < 8M to within our numerical error (< 10 -3 ). 



